#############################################
## PROGRAM NAME: 203e_res_msa_noimp.R      ##
## AUTHOR: MATT MLECZKO                    ##
## INPUTS:                                 ##
##         201_af_msa_fin.Rda              ##
##         202a_m1_msa_noimp.Rda -         ##
##         202a_m12_msa_noimp.Rda          ##
##                                         ##
## OUTPUTS:                                ##
##                                         ## 
## PURPOSE: Plot MSA no imp results        ##
##                                         ##
## LIST OF UPDATES:                        ##
#############################################

#log <- file(# USER DEFINED PATH AND FILE NAME HERE #)
#sink(log, append=TRUE)
#sink(log, append=TRUE, type="message")

## load libraries ##
library(tidyverse)

data_path <- # USER DEFINED PATH HERE #
setwd(data_path)

load("201_af_msa_fin.Rda")

load("202a_m1_msa_noimp.Rda")
load("202a_m2_msa_noimp.Rda")
load("202a_m3_msa_noimp.Rda")
load("202a_m4_msa_noimp.Rda")
load("202a_m5_msa_noimp.Rda")
load("202a_m6_msa_noimp.Rda")
load("202a_m7_msa_noimp.Rda")
load("202a_m8_msa_noimp.Rda")
load("202a_m9_msa_noimp.Rda")
load("202a_m10_msa_noimp.Rda")
load("202a_m11_msa_noimp.Rda")
load("202a_m12_msa_noimp.Rda")


msa.res1 <- msa.res.noimp1
msa.res2 <- msa.res.noimp2
msa.res3 <- msa.res.noimp3
msa.res4 <- msa.res.noimp4
msa.res5 <- msa.res.noimp5
msa.res6 <- msa.res.noimp6
msa.res7 <- msa.res.noimp7
msa.res8 <- msa.res.noimp8
msa.res9 <- msa.res.noimp9
msa.res10 <- msa.res.noimp10
msa.res11 <- msa.res.noimp11
msa.res12 <- msa.res.noimp12

## count non-missing obs ##

sum(!is.na(af.msa.fin$mgr_npov_2022))
sum(!is.na(af.msa.fin$mgr_hpov_2022))
sum(!is.na(af.msa.fin$mgr_hpov_cc_2022))

sum(!is.na(af.msa.fin$rb_npov_2022))
sum(!is.na(af.msa.fin$rb_hpov_2022))
sum(!is.na(af.msa.fin$rb_hpov_cc_2022))

sum(!is.na(af.msa.fin$mpv_npov_2022))
sum(!is.na(af.msa.fin$mpv_hpov_2022))
sum(!is.na(af.msa.fin$mpv_hpov_cc_2022))

sum(!is.na(af.msa.fin$eli_ushare_2022))
sum(!is.na(af.msa.fin$vli_ushare_2022))
sum(!is.na(af.msa.fin$li_ushare_2022))

## compile the results ## 

ex.msm.res <- function(i, inlist) {
  
  res <- inlist[[i]][1]
  t <- as.data.frame(res[[1]]$coefficients[2,1])
  names(t) <- "estimate"
  
  return(t)
  
}

ex.noadj.res <- function(i, inlist) {
  
  res <- inlist[[i]][3]
  t <- as.data.frame(res[[1]]$coefficients[2,1])
  names(t) <- "estimate"
  
  return(t)
  
}

ex.adl.res <- function(i, inlist) {
  
  res <- inlist[[i]][4]
  t <- as.data.frame(res[[1]]$coefficients[2,1])
  names(t) <- "estimate"
  
  return(t)
  
}


ex.fe.res <- function(i, inlist) {
  
  res <- inlist[[i]][5]
  t <- as.data.frame(res[[1]]$coefficients[2,1])
  names(t) <- "estimate"
  
  return(t)
  
}

## msm ##

m1.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res1)

m1.msm.res.ov <- bind_rows(m1.msm.res)


m2.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res2)

m2.msm.res.ov <- bind_rows(m2.msm.res)

m3.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res3)

m3.msm.res.ov <- bind_rows(m3.msm.res) 

m4.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res4)

m4.msm.res.ov <- bind_rows(m4.msm.res) 


m5.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res5)

m5.msm.res.ov <- bind_rows(m5.msm.res) 

m6.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res6)

m6.msm.res.ov <- bind_rows(m6.msm.res) 

m7.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res7)

m7.msm.res.ov <- bind_rows(m7.msm.res) 

m8.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res8)

m8.msm.res.ov <- bind_rows(m8.msm.res)

m9.msm.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = msa.res9)

m9.msm.res.ov <- bind_rows(m9.msm.res)

m10.msm.res <- lapply(seq(1,100,by=1),
                      ex.msm.res,
                      inlist = msa.res10)

m10.msm.res.ov <- bind_rows(m10.msm.res) 

m11.msm.res <- lapply(seq(1,100,by=1),
                      ex.msm.res,
                      inlist = msa.res11)

m11.msm.res.ov <- bind_rows(m11.msm.res)

m12.msm.res <- lapply(seq(1,100,by=1),
                      ex.msm.res,
                      inlist = msa.res12)

m12.msm.res.ov <- bind_rows(m12.msm.res)


## no adj ## 

m1.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res1)

m1.noadj.res.ov <- bind_rows(m1.noadj.res)


m2.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res2)

m2.noadj.res.ov <- bind_rows(m2.noadj.res)

m3.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res3)

m3.noadj.res.ov <- bind_rows(m3.noadj.res) 

m4.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res4)

m4.noadj.res.ov <- bind_rows(m4.noadj.res) 


m5.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res5)

m5.noadj.res.ov <- bind_rows(m5.noadj.res) 

m6.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res6)

m6.noadj.res.ov <- bind_rows(m6.noadj.res) 

m7.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res7)

m7.noadj.res.ov <- bind_rows(m7.noadj.res) 

m8.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res8)

m8.noadj.res.ov <- bind_rows(m8.noadj.res)

m9.noadj.res  <- lapply(seq(1,100,by=1),
                     ex.noadj.res,
                     inlist = msa.res9)

m9.noadj.res.ov <- bind_rows(m9.noadj.res)

m10.noadj.res  <- lapply(seq(1,100,by=1),
                      ex.noadj.res,
                      inlist = msa.res10)

m10.noadj.res.ov <- bind_rows(m10.noadj.res) 

m11.noadj.res  <- lapply(seq(1,100,by=1),
                      ex.noadj.res,
                      inlist = msa.res11)

m11.noadj.res.ov <- bind_rows(m11.noadj.res)

m12.noadj.res  <- lapply(seq(1,100,by=1),
                      ex.noadj.res,
                      inlist = msa.res12)

m12.noadj.res.ov <- bind_rows(m12.noadj.res)

## adl ## 

m1.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res1)

m1.adl.res.ov <- bind_rows(m1.adl.res)


m2.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res2)

m2.adl.res.ov <- bind_rows(m2.adl.res)

m3.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res3)

m3.adl.res.ov <- bind_rows(m3.adl.res) 

m4.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res4)

m4.adl.res.ov <- bind_rows(m4.adl.res) 


m5.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res5)

m5.adl.res.ov <- bind_rows(m5.adl.res) 

m6.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res6)

m6.adl.res.ov <- bind_rows(m6.adl.res) 

m7.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res7)

m7.adl.res.ov <- bind_rows(m7.adl.res) 

m8.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res8)

m8.adl.res.ov <- bind_rows(m8.adl.res)

m9.adl.res  <- lapply(seq(1,100,by=1),
                        ex.adl.res,
                        inlist = msa.res9)

m9.adl.res.ov <- bind_rows(m9.adl.res)

m10.adl.res  <- lapply(seq(1,100,by=1),
                         ex.adl.res,
                         inlist = msa.res10)

m10.adl.res.ov <- bind_rows(m10.adl.res) 

m11.adl.res  <- lapply(seq(1,100,by=1),
                         ex.adl.res,
                         inlist = msa.res11)

m11.adl.res.ov <- bind_rows(m11.adl.res)

m12.adl.res  <- lapply(seq(1,100,by=1),
                         ex.adl.res,
                         inlist = msa.res12)

m12.adl.res.ov <- bind_rows(m12.adl.res)


## fe ## 

m1.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res1)

m1.fe.res.ov <- bind_rows(m1.fe.res)


m2.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res2)

m2.fe.res.ov <- bind_rows(m2.fe.res)

m3.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res3)

m3.fe.res.ov <- bind_rows(m3.fe.res) 

m4.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res4)

m4.fe.res.ov <- bind_rows(m4.fe.res) 


m5.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res5)

m5.fe.res.ov <- bind_rows(m5.fe.res) 

m6.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res6)

m6.fe.res.ov <- bind_rows(m6.fe.res) 

m7.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res7)

m7.fe.res.ov <- bind_rows(m7.fe.res) 

m8.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res8)

m8.fe.res.ov <- bind_rows(m8.fe.res)

m9.fe.res  <- lapply(seq(1,100,by=1),
                      ex.fe.res,
                      inlist = msa.res9)

m9.fe.res.ov <- bind_rows(m9.fe.res)

m10.fe.res  <- lapply(seq(1,100,by=1),
                       ex.fe.res,
                       inlist = msa.res10)

m10.fe.res.ov <- bind_rows(m10.fe.res) 

m11.fe.res  <- lapply(seq(1,100,by=1),
                       ex.fe.res,
                       inlist = msa.res11)

m11.fe.res.ov <- bind_rows(m11.fe.res)

m12.fe.res  <- lapply(seq(1,100,by=1),
                       ex.fe.res,
                       inlist = msa.res12)

m12.fe.res.ov <- bind_rows(m12.fe.res)


## plot the results ##

msm.res.in <- data.frame(var = c("Median gross rent (nat)",
                                 "Median gross rent (hpov)",
                                 "Median gross rent (hpov cc)",
                                 "Rent burden (nat)",
                                 "Rent burden (hpov)",
                                 "Rent burden (hpov cc)",
                                 "MPV (nat)",
                                 "MPV (hpov)",
                                 "MPV (hpov cc)",
                                 "ELI supply",
                                 "VLI supply",
                                 "LI supply"),
                         coeff = c(mean(m1.msm.res.ov$estimate),
                                   mean(m2.msm.res.ov$estimate),
                                   mean(m3.msm.res.ov$estimate),
                                   mean(m4.msm.res.ov$estimate),
                                   mean(m5.msm.res.ov$estimate),
                                   mean(m6.msm.res.ov$estimate),
                                   mean(m7.msm.res.ov$estimate),
                                   mean(m8.msm.res.ov$estimate),
                                   mean(m9.msm.res.ov$estimate),
                                   mean(m10.msm.res.ov$estimate),
                                   mean(m11.msm.res.ov$estimate),
                                   mean(m12.msm.res.ov$estimate)),
                         lb = c(quantile(m1.msm.res.ov$estimate,0.025),
                                quantile(m2.msm.res.ov$estimate,0.025),
                                quantile(m3.msm.res.ov$estimate,0.025),
                                quantile(m4.msm.res.ov$estimate,0.025),
                                quantile(m5.msm.res.ov$estimate,0.025),
                                quantile(m6.msm.res.ov$estimate,0.025),
                                quantile(m7.msm.res.ov$estimate,0.025),
                                quantile(m8.msm.res.ov$estimate,0.025),
                                quantile(m9.msm.res.ov$estimate,0.025),
                                quantile(m10.msm.res.ov$estimate,0.025),
                                quantile(m11.msm.res.ov$estimate,0.025),
                                quantile(m12.msm.res.ov$estimate,0.025)),
                         ub = c(quantile(m1.msm.res.ov$estimate,0.975),
                                quantile(m2.msm.res.ov$estimate,0.975),
                                quantile(m3.msm.res.ov$estimate,0.975),
                                quantile(m4.msm.res.ov$estimate,0.975),
                                quantile(m5.msm.res.ov$estimate,0.975),
                                quantile(m6.msm.res.ov$estimate,0.975),
                                quantile(m7.msm.res.ov$estimate,0.975),
                                quantile(m8.msm.res.ov$estimate,0.975),
                                quantile(m9.msm.res.ov$estimate,0.975),
                                quantile(m10.msm.res.ov$estimate,0.975),
                                quantile(m11.msm.res.ov$estimate,0.975),
                                quantile(m12.msm.res.ov$estimate,0.975)))

## no adj ## 

noadj.res.in <- data.frame(var = c("Median gross rent (nat)",
                                   "Median gross rent (hpov)",
                                   "Median gross rent (hpov cc)",
                                   "Rent burden (nat)",
                                   "Rent burden (hpov)",
                                   "Rent burden (hpov cc)",
                                   "MPV (nat)",
                                   "MPV (hpov)",
                                   "MPV (hpov cc)",
                                   "ELI supply",
                                   "VLI supply",
                                   "LI supply"),
                           coeff = c(mean(m1.noadj.res.ov$estimate),
                                     mean(m2.noadj.res.ov$estimate),
                                     mean(m3.noadj.res.ov$estimate),
                                     mean(m4.noadj.res.ov$estimate),
                                     mean(m5.noadj.res.ov$estimate),
                                     mean(m6.noadj.res.ov$estimate),
                                     mean(m7.noadj.res.ov$estimate),
                                     mean(m8.noadj.res.ov$estimate),
                                     mean(m9.noadj.res.ov$estimate),
                                     mean(m10.noadj.res.ov$estimate),
                                     mean(m11.noadj.res.ov$estimate),
                                     mean(m12.noadj.res.ov$estimate)),
                           lb = c(quantile(m1.noadj.res.ov$estimate,0.025),
                                  quantile(m2.noadj.res.ov$estimate,0.025),
                                  quantile(m3.noadj.res.ov$estimate,0.025),
                                  quantile(m4.noadj.res.ov$estimate,0.025),
                                  quantile(m5.noadj.res.ov$estimate,0.025),
                                  quantile(m6.noadj.res.ov$estimate,0.025),
                                  quantile(m7.noadj.res.ov$estimate,0.025),
                                  quantile(m8.noadj.res.ov$estimate,0.025),
                                  quantile(m9.noadj.res.ov$estimate,0.025),
                                  quantile(m10.noadj.res.ov$estimate,0.025),
                                  quantile(m11.noadj.res.ov$estimate,0.025),
                                  quantile(m12.noadj.res.ov$estimate,0.025)),
                           ub = c(quantile(m1.noadj.res.ov$estimate,0.975),
                                  quantile(m2.noadj.res.ov$estimate,0.975),
                                  quantile(m3.noadj.res.ov$estimate,0.975),
                                  quantile(m4.noadj.res.ov$estimate,0.975),
                                  quantile(m5.noadj.res.ov$estimate,0.975),
                                  quantile(m6.noadj.res.ov$estimate,0.975),
                                  quantile(m7.noadj.res.ov$estimate,0.975),
                                  quantile(m8.noadj.res.ov$estimate,0.975),
                                  quantile(m9.noadj.res.ov$estimate,0.975),
                                  quantile(m10.noadj.res.ov$estimate,0.975),
                                  quantile(m11.noadj.res.ov$estimate,0.975),
                                  quantile(m12.noadj.res.ov$estimate,0.975)))

## adl ## 

adl.res.in <- data.frame(var = c("Median gross rent (nat)",
                                 "Median gross rent (hpov)",
                                 "Median gross rent (hpov cc)",
                                 "Rent burden (nat)",
                                 "Rent burden (hpov)",
                                 "Rent burden (hpov cc)",
                                 "MPV (nat)",
                                 "MPV (hpov)",
                                 "MPV (hpov cc)",
                                 "ELI supply",
                                 "VLI supply",
                                 "LI supply"),
                         coeff = c(mean(m1.adl.res.ov$estimate),
                                   mean(m2.adl.res.ov$estimate),
                                   mean(m3.adl.res.ov$estimate),
                                   mean(m4.adl.res.ov$estimate),
                                   mean(m5.adl.res.ov$estimate),
                                   mean(m6.adl.res.ov$estimate),
                                   mean(m7.adl.res.ov$estimate),
                                   mean(m8.adl.res.ov$estimate),
                                   mean(m9.adl.res.ov$estimate),
                                   mean(m10.adl.res.ov$estimate),
                                   mean(m11.adl.res.ov$estimate),
                                   mean(m12.adl.res.ov$estimate)),
                         lb = c(quantile(m1.adl.res.ov$estimate,0.025),
                                quantile(m2.adl.res.ov$estimate,0.025),
                                quantile(m3.adl.res.ov$estimate,0.025),
                                quantile(m4.adl.res.ov$estimate,0.025),
                                quantile(m5.adl.res.ov$estimate,0.025),
                                quantile(m6.adl.res.ov$estimate,0.025),
                                quantile(m7.adl.res.ov$estimate,0.025),
                                quantile(m8.adl.res.ov$estimate,0.025),
                                quantile(m9.adl.res.ov$estimate,0.025),
                                quantile(m10.adl.res.ov$estimate,0.025),
                                quantile(m11.adl.res.ov$estimate,0.025),
                                quantile(m12.adl.res.ov$estimate,0.025)),
                         ub = c(quantile(m1.adl.res.ov$estimate,0.975),
                                quantile(m2.adl.res.ov$estimate,0.975),
                                quantile(m3.adl.res.ov$estimate,0.975),
                                quantile(m4.adl.res.ov$estimate,0.975),
                                quantile(m5.adl.res.ov$estimate,0.975),
                                quantile(m6.adl.res.ov$estimate,0.975),
                                quantile(m7.adl.res.ov$estimate,0.975),
                                quantile(m8.adl.res.ov$estimate,0.975),
                                quantile(m9.adl.res.ov$estimate,0.975),
                                quantile(m10.adl.res.ov$estimate,0.975),
                                quantile(m11.adl.res.ov$estimate,0.975),
                                quantile(m12.adl.res.ov$estimate,0.975)))


fe.res.in <- data.frame(var = c("Median gross rent (nat)",
                                "Median gross rent (hpov)",
                                "Median gross rent (hpov cc)",
                                "Rent burden (nat)",
                                "Rent burden (hpov)",
                                "Rent burden (hpov cc)",
                                "MPV (nat)",
                                "MPV (hpov)",
                                "MPV (hpov cc)",
                                "ELI supply",
                                "VLI supply",
                                "LI supply"),
                        coeff = c(mean(m1.fe.res.ov$estimate),
                                  mean(m2.fe.res.ov$estimate),
                                  mean(m3.fe.res.ov$estimate),
                                  mean(m4.fe.res.ov$estimate),
                                  mean(m5.fe.res.ov$estimate),
                                  mean(m6.fe.res.ov$estimate),
                                  mean(m7.fe.res.ov$estimate),
                                  mean(m8.fe.res.ov$estimate),
                                  mean(m9.fe.res.ov$estimate),
                                  mean(m10.fe.res.ov$estimate),
                                  mean(m11.fe.res.ov$estimate),
                                  mean(m12.fe.res.ov$estimate)),
                        lb = c(quantile(m1.fe.res.ov$estimate,0.025),
                               quantile(m2.fe.res.ov$estimate,0.025),
                               quantile(m3.fe.res.ov$estimate,0.025),
                               quantile(m4.fe.res.ov$estimate,0.025),
                               quantile(m5.fe.res.ov$estimate,0.025),
                               quantile(m6.fe.res.ov$estimate,0.025),
                               quantile(m7.fe.res.ov$estimate,0.025),
                               quantile(m8.fe.res.ov$estimate,0.025),
                               quantile(m9.fe.res.ov$estimate,0.025),
                               quantile(m10.fe.res.ov$estimate,0.025),
                               quantile(m11.fe.res.ov$estimate,0.025),
                               quantile(m12.fe.res.ov$estimate,0.025)),
                        ub = c(quantile(m1.fe.res.ov$estimate,0.975),
                               quantile(m2.fe.res.ov$estimate,0.975),
                               quantile(m3.fe.res.ov$estimate,0.975),
                               quantile(m4.fe.res.ov$estimate,0.975),
                               quantile(m5.fe.res.ov$estimate,0.975),
                               quantile(m6.fe.res.ov$estimate,0.975),
                               quantile(m7.fe.res.ov$estimate,0.975),
                               quantile(m8.fe.res.ov$estimate,0.975),
                               quantile(m9.fe.res.ov$estimate,0.975),
                               quantile(m10.fe.res.ov$estimate,0.975),
                               quantile(m11.fe.res.ov$estimate,0.975),
                               quantile(m12.fe.res.ov$estimate,0.975)))



## set 1: median rents ## 

msm.fin.res1 <- msm.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above-average (N=349)",
                           var == "Median gross rent (hpov)" ~ "High-poverty (N=221)",
                           var == "Median gross rent (hpov cc)" ~ "High-poverty cc (N=207)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 1 (N=349)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 1 (N=221)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 1 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res1 <- noadj.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above-average (N=349)",
                           var == "Median gross rent (hpov)" ~ "High-poverty (N=221)",
                           var == "Median gross rent (hpov cc)" ~ "High-poverty cc (N=207)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 2 (N=349)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 2 (N=221)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 2 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res1 <- adl.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above-average (N=349)",
                           var == "Median gross rent (hpov)" ~ "High-poverty (N=221)",
                           var == "Median gross rent (hpov cc)" ~ "High-poverty cc (N=207)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 3 (N=349)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 3 (N=221)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 3 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "ADL")

fe.fin.res1 <- fe.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above-average (N=349)",
                           var == "Median gross rent (hpov)" ~ "High-poverty (N=221)",
                           var == "Median gross rent (hpov cc)" ~ "High-poverty cc (N=207)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 4 (N=349)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 4 (N=221)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 4 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res1 <- rbind(msm.fin.res1,
                  noadj.fin.res1,
                  adl.fin.res1,
                  fe.fin.res1)

fin.res1$Model <- factor(fin.res1$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))


ggplot(fin.res1, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  theme_bw(base_size = 14) + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  theme(panel.spacing = unit(0, "lines"), 
        strip.background = element_blank(),
        strip.placement = "outside") +
  ggtitle("") + 
  xlab("") + 
  ylab("Estimate ($)") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))


## set 2: rent burden ##

msm.fin.res2 <- msm.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above-average (N=349)",
                           var == "Rent burden (hpov)" ~ "High-poverty (N=221)",
                           var == "Rent burden (hpov cc)" ~ "High-poverty cc (N=207)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 1 (N=349)",
                         var == "Rent burden (hpov)" ~ "High-poverty 1 (N=221)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 1 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res2 <- noadj.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above-average (N=349)",
                           var == "Rent burden (hpov)" ~ "High-poverty (N=221)",
                           var == "Rent burden (hpov cc)" ~ "High-poverty cc (N=207)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 2 (N=349)",
                         var == "Rent burden (hpov)" ~ "High-poverty 2 (N=221)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 2 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res2 <- adl.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above-average (N=349)",
                           var == "Rent burden (hpov)" ~ "High-poverty (N=221)",
                           var == "Rent burden (hpov cc)" ~ "High-poverty cc (N=207)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 3 (N=349)",
                         var == "Rent burden (hpov)" ~ "High-poverty 3 (N=221)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 3 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "ADL")

fe.fin.res2 <- fe.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above-average (N=349)",
                           var == "Rent burden (hpov)" ~ "High-poverty (N=221)",
                           var == "Rent burden (hpov cc)" ~ "High-poverty cc (N=207)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 4 (N=349)",
                         var == "Rent burden (hpov)" ~ "High-poverty 4 (N=221)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 4 (N=207)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res2 <- rbind(msm.fin.res2,
                  noadj.fin.res2,
                  adl.fin.res2,
                  fe.fin.res2)

fin.res2$Model <- factor(fin.res2$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))


ggplot(fin.res2, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  theme_bw(base_size = 14) + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  theme(panel.spacing = unit(0, "lines"), 
        strip.background = element_blank(),
        strip.placement = "outside") +
  ggtitle("") + 
  xlab("") + 
  ylab("Estimate") +
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))

## set 3: MPV in impoverished neighborhoods ##

msm.fin.res3 <- msm.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above-average (N=349)",
                           var == "MPV (hpov)" ~ "High-poverty (N=218)",
                           var == "MPV (hpov cc)" ~ "High-poverty cc (N=204)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 1 (N=349)",
                         var == "MPV (hpov)" ~ "High-poverty 1 (N=218)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 1 (N=204)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res3 <- noadj.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above-average (N=349)",
                           var == "MPV (hpov)" ~ "High-poverty (N=218)",
                           var == "MPV (hpov cc)" ~ "High-poverty cc (N=204)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 2 (N=349)",
                         var == "MPV (hpov)" ~ "High-poverty 2 (N=218)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 2 (N=204)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res3 <- adl.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above-average (N=349)",
                           var == "MPV (hpov)" ~ "High-poverty (N=218)",
                           var == "MPV (hpov cc)" ~ "High-poverty cc (N=204)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 3 (N=349)",
                         var == "MPV (hpov)" ~ "High-poverty 3 (N=218)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 3 (N=204)"),
         lb = lb,
         ub = ub,
         Model = "ADL")


fe.fin.res3 <- fe.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above-average (N=349)",
                           var == "MPV (hpov)" ~ "High-poverty (N=218)",
                           var == "MPV (hpov cc)" ~ "High-poverty cc (N=204)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 4 (N=349)",
                         var == "MPV (hpov)" ~ "High-poverty 4 (N=218)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 4 (N=204)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.res3 <- rbind(msm.fin.res3,
                  noadj.fin.res3,
                  adl.fin.res3,
                  fe.fin.res3)

fin.res3$Model <- factor(fin.res3$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))

ggplot(fin.res3, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  theme_bw(base_size = 14) + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  theme(panel.spacing = unit(0, "lines"), 
        strip.background = element_blank(),
        strip.placement = "outside") +
  ggtitle("") + 
  xlab("") + 
  ylab("Estimate ($)") +
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))

## set 4: Supply ##

msm.fin.res4 <- msm.res.in %>%
  filter(var %in% c("ELI supply",
                    "VLI supply",
                    "LI supply")) %>%
  mutate(group = case_when(var == "ELI supply" ~ "ELI supply (N=349)",
                           var == "VLI supply" ~ "VLI supply (N=349)",
                           var == "LI supply" ~ "LI supply (N=349)"),
         var = case_when(var == "ELI supply" ~ "ELI supply 1 (N=349)",
                         var == "VLI supply" ~ "VLI supply 1 (N=349)",
                         var == "LI supply" ~ "LI supply 1 (N=349)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.res4 <- noadj.res.in %>%
  filter(var %in% c("ELI supply",
                    "VLI supply",
                    "LI supply")) %>%
  mutate(group = case_when(var == "ELI supply" ~ "ELI supply (N=349)",
                           var == "VLI supply" ~ "VLI supply (N=349)",
                           var == "LI supply" ~ "LI supply (N=349)"),
         var = case_when(var == "ELI supply" ~ "ELI supply 2 (N=349)",
                         var == "VLI supply" ~ "VLI supply 2 (N=349)",
                         var == "LI supply" ~ "LI supply 2 (N=349)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.res4 <- adl.res.in %>%
  filter(var %in% c("ELI supply",
                    "VLI supply",
                    "LI supply")) %>%
  mutate(group = case_when(var == "ELI supply" ~ "ELI supply (N=349)",
                           var == "VLI supply" ~ "VLI supply (N=349)",
                           var == "LI supply" ~ "LI supply (N=349)"),
         var = case_when(var == "ELI supply" ~ "ELI supply 3 (N=349)",
                         var == "VLI supply" ~ "VLI supply 3 (N=349)",
                         var == "LI supply" ~ "LI supply 3 (N=349)"),
         lb = lb,
         ub = ub,
         Model = "ADL")


fe.fin.res4 <- fe.res.in %>%
  filter(var %in% c("ELI supply",
                    "VLI supply",
                    "LI supply")) %>%
  mutate(group = case_when(var == "ELI supply" ~ "ELI supply (N=349)",
                           var == "VLI supply" ~ "VLI supply (N=349)",
                           var == "LI supply" ~ "LI supply (N=349)"),
         var = case_when(var == "ELI supply" ~ "ELI supply 4 (N=349)",
                         var == "VLI supply" ~ "VLI supply 4 (N=349)",
                         var == "LI supply" ~ "LI supply 4 (N=349)"),
         lb = lb,
         ub = ub,
         Model = "FE")


fin.res4 <- rbind(msm.fin.res4,
                  noadj.fin.res4,
                  adl.fin.res4,
                  fe.fin.res4)

fin.res4$group <- factor(fin.res4$group,
                         levels = c("ELI supply (N=349)",
                                    "VLI supply (N=349)",
                                    "LI supply (N=349)"))

fin.res4$Model <- factor(fin.res4$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))

ggplot(fin.res4, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  theme_bw(base_size = 14) + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  theme(panel.spacing = unit(0, "lines"), 
        strip.background = element_blank(),
        strip.placement = "outside") +
  ggtitle("") + 
  xlab("") + 
  ylab("Estimate") +
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))


## END OF PROGRAM ##

#sink()
